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We propose a theoretical framework for reconstructing tip-surface interactions using the inter- 
modulation technique when more than one eigenmode is required to describe the cantilever motion. 
Two particular cases of bimodal motion are studied numerically: one bending and one torsional 
mode, and two bending modes. We demonstrate the possibility of accurate reconstruction of a 
two-dimensional conservative force field for the former case, while dissipative forces are studied for 
the latter. 
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I. INTRODUCTION 

Atomic force microscopy^ (AFM) has become one of 
the most important tools for the study of nanometer- 
scale surface properties of a wide range of materials. 
The initial goal of AFM was surface topography imag- 
ing which was performed by scanning a cantilever with 
a sharp tip over a surface while keeping its deflection 
constant. It was later realized that the reconstruction of 
tip-surface interactions was possible and that these in- 
teractions contain valuable information about material 
properties^. One of the first reconstruction methods 
was based on measurement of the quasi-static bending 
of the cantilever beam as its base was slowly moved to- 
ward and away from the surface. Two drawbacks of this 
method are the slow speed of measurement and the lack 
of ability to reconstruct dissipative forces which are al- 
ways present in tip-surface interactions due to non-elastic 
deformations of the sample, breaking chemical bonds, or 
other irreversible processes^^. 

The development of dynamic AFM opened new path- 
ways for a more profound study of tip-surface interac- 
tions. In dynamic AFM the cantilever is treated as 
an undcrdamped oscillator (high Q-factor) driven at a 
resonance where the response to external forces is en- 
hanced by a factor Q. A large number of methods 
for determining the tip-surface interaction have been 
devisecP^, some making use of amplitude or frequency 
modulatio n 16 " 25 !. In this paper, we consider intermodu- 
lation AFM (ImAFM^H which is unique in its ability to 
rapidly extract a large amount of information abou t tip- 
surface interactions in a simple and convenient wajE^U. 
With ImAFM we have the possibility of high resolution 
surface property mapping at interactive scan speeds and 
reconstruction of the tip-surface interaction at each pixel. 



The main idea underlying ImAFM is to use the non- 
linear tip-surface forces to create high-order intermod- 
ulation of discrete tones in a frequency comb®. The 
method can be generally applied to any resonator sub- 
ject to a nonlinear force when it is driven with at least 
two frequencies u>di and u>d2- The nonlinear response 
in frequency domain occurs not only at the drive fre- 
quencies and their harmonics nujdi and muid2, where n 
and m are integers, but also at their linear combinations 
u) = nujdi + mujd2 or intermodulation products. Due to 
the signal enhancement near resonance and finite detec- 
tor noise of the measurement, a typical spectrum of can- 
tilever motion can be obtained only in narrow frequency 
band near resonance. Concentrating as many intermod- 
ulation products as possible in this resonant detection 
band results in more information for reconstruction of 
the nonlinear forces (fig. [TJ). 

The sensitivity enhancement of dynamic AFM occurs 
not only at one single eigenmode, but also at several 
eigenmodes simultaneously. Measurement of the sam- 
ple response at higher frequencies allows for better ma- 
terial imaging contrasts and improves accuracy of the 
tip-surface force reconstructiorl22H2lD if we a i so excite a 
torsional resonance of the cantilever we open up the pos- 
sibility of measuring lateral forces acting on the tip^HSSl 
Additionally, when driving the cantilever at more than 
one eigenmode, new frequency bands become available 
for collecting intermodulation products, resulting from 
the nonlinear force which couples multiple eigenmodes^. 
While ImAF M has b een well studied for the case of one 
flexural mode^^iQIIHlS ^he mu lti-eigenmode problem re- 
mains open. In the current investigation we explore the 
possibility of reconstructing two-dimensional tip-surface 
force fields using two collinear eigenmodes (e.g. two flex- 
ural or two torsional modes) as well as two orthogonal 
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FIG. 1. (Color online) Appearance of new frequencies in the 
spectrum of the tip motion in ImAFM due to the nonlinear- 
ities in tip-surface interaction. The dashed line shows the 
characteristic range of the tip-surface force. 



modes (e.g. one flexural and one torsional mode). 

The paper is organized as follows: in section [II] we 
consider general properties of the cantilever model and 
tip-surface interactions. In section |III| we recount the 
main principles of ImAFM, demonstrate how to obtain 
intermodulation spectra of tip-surface forces , and develop 
an extension of the spectral fitting methooP^SI of force 
reconstruction for the multimodal case. In Section llVl nu- 
merical results for the reconstruction of two-dimensional 
conservative and dissipative forces are presented. Section 
|V| concludes with a discussion and summary. 



II. MODEL 

In order to proceed with the multimodal problem 
we should start from some general discussion of can- 
tilever dynamic j 44 l 45 | If we are interested in study- 
ing only flexural modes, one-dimensional Euler-Bernoulli 
beam theorjSsHU] is sufficient. Incorporating torsional 
modes requires the cantilever to be regarded as a two- 
dimensional object. Many theories have been developed 
describing the continuum mechanics of two-dimensional 
plategS2HS3|. A general description of an arbitrary two- 
dimensional cantilever (fig. [2| is given by the governing 
equation 



{Q xy + Qt)[w (x,y,t)]=F + i 



(1) 



with an appropriate set of boundary conditions. Here x 
and y are the two space coordinates, t is time, w{x,y,i) 
is a deflection normal to the x — y plane of the plate at 
rest, Gxy is a space coupling operator which corresponds 
to the two-dimensional mathematical model used to de- 
termine the stresses and deformations in thin plates, and 
Qt is a time evolution operator which represents inertia 
and damping. F and f are a scalar quantities which are 



projections on to w, of the two-dimensional vector force 
field acting on the tip, and the drive, respectively. In 
general, the force F can depend explicitly on the deflec- 
tion w, its velocity w, past trajectories {w,w}\_ 00 and 
time t. In this paper we restrict ourselves to F which 
are not dependent on past trajectories. We consider only 
the case of small deflections, so Q xy and Gt are linear. 
For example, within the Kirchhoff-Love plate theory^ 
the space operator for a homogeneous cantilever reads 



2h 3 E 
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(2) 



where 2h is the cantilever thickness, E the Young's mod- 
ulus and v the Poisson's ratio. The theory assumes that 
a mid-surface plane can be used to represent a three- 
dimensional plate in a two-dimensional form. The time 
evolution operator Q t consists of an inertial term, and 
for the case of a homogeneous viscous media, a linear 
damping term 



Qt ■= 2ph 



dt 2 



•27 



d_ 
dt 



(3) 



where p is the cantilever density and 7 a damping coeffi- 
cient. 

The linearity of Q xy and Q t gives a dynamics of the two- 
dimensional cantilever which is well approximated by a 
system of differential equations for the generalized coor- 
dinates qi representing deflections of its normal modes 
(see appendix | for a derivation) 
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F t (t) + U (t) (4) 



Here each generalized coordinate qi has a stiffness fcj , res- 
onant frequency u>i and quality factor Qi . F^ (t) is the pro- 
jection of the nonlinear force and fj(i) is the projection 
of the drive force on to the i th eigenmode. The problem 
of mapping the eigenmode coordinates qi onto the phys- 
ical position of the tip r = (z,y) J , (hereafter j denotes 
transpose) as well as the relationship between Fi and the 
actual vector force field acting on the tip F ts (r, r), is de- 
termined by the cantilever and tip geometry as briefly 
discussed in section ITVl 



III. RECONSTRUCTING FORCE FROM 
INTERMODULATION WITH MULTIPLE 
EIGENMODES 

We recall the basic steps to get information about Fi 
in the frequency domain using the intermodulation tech- 
nique. Firstly, we excite one or several generalized co- 
ordinates Q with the drive fj and measure its motion 
spectrum at the height ho, far from the sample surface, 
so Fj =0 for all t. This free oscillation spectrum is 

ho 

denoted 



qi\ ho = Xdi 



(5) 
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FIG. 2. (Color online) Schematic representation of the two-dimensional cantilever and detection principle in AFM. 



where the hat represents the discrete Fourier transform 
(DFT) 

OO 

g(«)=^[«(*)]:= E <l(t)e- 2mm5ut (6) 

m— — oo 

and the linear response function \i is introduced 



Xi (uj) = k t 1 



Qi V w 



(7) 



Using pairs of drive tones with the same amplitude A t 
which are separated by Slo and placed close to each res- 
onance riibuj < LJi < (rii + 1)Suj, where rii is an integer, 

f t (i) = Ai [cos (mSuit) + cos ((rii + I) Suit)] (8) 



we obtain a free, linear response spectrum (eq. ([5])) con- 
sisting components only at the drive frequencies Wj±&j/2 
(fig.[TJ). Data acquisition in ImAFM should be performed 
over at least one total period of the drive T = 2it/5uj 
while measuring subsequent perio ds only improves the 
SNR of the measured spectre^^ES. 

We then we move the cantilever closer to the sample to 
the engaged height h, so the oscillating tip starts to feel 
interaction with a surface (but not too close, so that there 
is zero static deflection of the cantilever) and measure the 
motion spectrum again 



nh 



X, 



(Pi + u) 



(9) 



Finally, the difference between this engaged motion spec- 
trum ([9]) and free motion spectrum ([5]) yields the desired 
interaction force spectrum 



i\h 



(10) 



where finite difference operator A is used for short. Thus, 
we have obtained information about Fi in the frequency 
space representation of the motion and we are in a po- 
sition to discuss the reconstruction of its dependence on 
generalized coordinates qi and velocities qi. 



Under ideal conditions, given the full difference spec- 
trum Acji and corresponding response function Xii it is 
possible to find Fj(t) as the inverse Fourier transform 
of Fi(oj) and then trivially recover its coordinate depen- 
dence Fi({qi, qi}) using the measured motion qi(t). How- 
ever, in real experiments this naive approach fails due to 
the strong frequency dependence of xi an d the measure- 
ment limitations imposed by detection noise. In prac- 
tice, almost all of the spectrum qi is buried under detec- 
tor noise except for a narrow band near its resonant fre- 
quency where the signal-to-noise ratio (SNR) meets the 
thermal limit. Usually the number of resolvable spectral 
components Bi in a band surrounding each eigenmode 
resonance is limited to a few dozen, depending on the 
difference frequency Su) and the forces experienced dur- 
ing the interaction. The use of several (iV) eigenmodes 
allows for a lager total number of frequency components 
B = J2iLi Bi for force reconstruction. 

Different methods have been elaborated for the single 
mode reconstruction problem using this limit ed amou nt 
of response in the resonant detection bancP^USHm 
Here we deve lop an extension of the spectral fitting 
methocP 7 * 28 * 43 ! for the multimodal problem as it allows 
for a straightforward generalization without involving 
any sophisticated concepts^. Following the method's 
main idea, one assumes a tip-surface force in the form of 
some known model function F^qi, . . . , q^; g) with a vec- 
tor g = (<?!, . . . gp) 1 of P unknown parameters. Fitting 

the calculated spectrum Fi{ui) to the measured Fi(u>) (eq. 
( |10[ )), we minimize the error function in the frequency 
domain, in a least square sense 



mm ei 
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Fi - P { Ql 



■ j In; 



(11) 



The model Fi can be a particular phenomenological ex- 
pression, for example the van der Waals-Derjagin-Muller- 
Toporov (vdW-DMT) forctP or its modificationEPMl. 
However, in the general case we do not know the exact 
form of the interaction and should choose some generic 
function structure, for instance a truncated Taylor ex- 
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pansion in the following polynomial form 

Pi Pn 



Fi(qi, 



,Qn) = ^2 "' 9ii...i N qi ■ ■ ■ q l N = q T s 



ii=0 



jjv=0 



(12) 

Here Pi is the degree of the polynomial in qi 7 q and g are 
vectors of basis functions and parameters respectively, 
each being of size P — JI^li Pi 

q=(l, 1i, 1\, 92, 9i 92, 9?92, •■•) T 



(5o 



510... i 520. 



5oi..., 5n..., 521. 



(13) 

Although the polynomial model is more universal, it 
usually contains a much larger number of unknown 
parameters^ which do not directly correspond to physi- 
cal properties of the material or surface. 

Inserting (12) into equation (10), we obtain a system 



of linear equations for the polynomial coefficients gi l ...i N 
which is conveniently represented in matrix notation 



g = H+F 



(14) 



where H is a B x P matrix with rows Hk = -?fe[q T ] (k 
is the corresponding component of the discrete Fourier 
transform of q), H + its Moore-Penrose pseudoinverse, 
and Fj a vector of size Note that total number of 

unknown parameters P can not be greater than B and 
in the case of an overdetermined system (P < B) the 
pseudoinverse will give a unique solution to 
least square sense (11) 
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Reconstruction of the ve 
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dependent part of the force F{ can be performed in the 
same way as for a conservative force, by polynomial ex- 
pansion in qi with corresponding coefficients. 



IV. RECONSTRUCTION OF 
TWO-DIMENSIONAL TIP-SURFACE FORCES 
FROM INTERMODULATION AFM SPECTRA 

Thus the reconstruction of tip-surface force from multi- 
ple eigenmodes is a straight-forward generalization of the 
single eigenmode problem, albeit with the complication 
of keeping track of multiple modes and the possibility 
that tones can intermodulate between these modes. From 
an algebraic point of view the spectral fitting method 
can be regarded as a multivariate interpolation 62 and the 
simplest model which is linea r in the parameters suffers 
from Runge's phenomena^^when high-order nonlinear- 
ities couple the multiple eigenmode coordinates. Further- 
more, reconstruction from many eigenmodes is a multi- 
dimensional problem with many model parameters and 
it will ultimately suffer from the need to either calibrate 
or determine these parameters from the limited number 
of intermodulation products that can be extracted from 
the narrow frequency bands near the resonances. 

In a real experiment, the AFM detector is typically 
not able measure over a frequency range which covers 
many eigenmodes of the cantilever. Furthermore, the 



AFM detector is only capable of measuring two signals, 
corresponding to two orthogonal motions of the can- 
tilever, that of flexing and twisting. The case of only two 
eigenmodes is therefore a reasonable simplification of the 
multi-modal problem which of great practical interest. In 
the following we restrict ourselves to this bimodal case. 

We are interested in reconstruction of the two- 
dimensional tip-surface vector force field F ts (r, r) which 
depends on the physical tip position in the y — z plane, 
r = (z, y) J . Before we proceed with the two-mode analy- 
sis we should map the set of forces Fj acting on the can- 
tilever eigenmodes qi onto the physical force F ts . With 
this aim, we can transform the basis set for defining qi to 
separate "pure" modes from "mixed" modes. Applying 
some transformation, the exact form of which depends 
on the geometrical shape of the cantilever, we obtain the 
pure eigencoordinates Zi and yi contributing only to the 
tip position perpendicular and parallel to the surface re- 
spectively. Here the two different Zi (or j/j) are collinear 
while Zi and yi are orthogonal. The remaining q* are 
mixed eigencoordinates of the cantilever, or cross-modes 
contributing to both coordinates z and y simultaneously, 
so that 



z = E z i + E q* • z 
v = E v% + E qf • y 



(15) 



Here the second terms are projections of the mixed eigen- 
coordinates onto the tip coordinate system (z,y). In do- 
ing so, the force F ts is projected onto the (z, y) so that its 
components parallel to the surface F z and perpendicular 
to the surface F y act on the corresponding pure modes, 
i.e. Fi = F z for each Zi, F{ — F y for each yi, and different 
force projections act on the mixed modes. 

Simultaneous excitation of the several pure eigenmodes 
coupled by nonlinear tip-surface interaction leads to exci- 
tation of the mixed modes which allows for measurements 
of the response in additional frequency bands. Although 
it would provide additional information, in this paper we 
investigate the simplest multimodal motion, that of pure 
bimodal motion without considering the cross-modes. In 
this case, analysis of the cantilever dynamics reduces to 
study of two characteristic regimes: bimodal motion of 
the collinear eigencoordinates (e.g. two flexural or tor- 
sional modes) and orthogonal eigencoordinates (one flex- 
ural and one torsional mode). Let us proceed with the 
former. 



A. Case 1: Two collinear modes 

This case corresponds to the dynamics of two flexu- 
ral modes z\ and z 2 , so that the total perpendicular tip 
deflection z = z\ + z 2 



fc i + *i) =F z (z,z)+h(t) 



;Z 2 



-z 2 + z 2 , 



F z (z,z)+t 2 (t) 



(16) 
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FIG. 3. (Color online) Spectrum of the engaged cantilever mo- 
tion with two fiexural modes. Two highest peaks are clearly 
seen above noise level near resonant frequencies ui and UJ2 
which consist of the intermodulation products depicted in in- 
sertions. The spectrum is obtained by integrating the system 
with the tip-surface force \22\ and consequent addition 
of a white noise. 



Here, on the right hand side we have the same vector 
component F z of the tip-surface force field depending on 
z and i. 

Firstly, let us consider a model for the position- 
dependent part of F z in some general form 



p z p, 
£(*i,*2) = EE 

i=0 j=0 



9ij z l z 2 



(17) 



This model requires determination of a large number 
PziPz — l)/2 of coefficients gij in order to define the 
polynomial of order P z = P Zl P Z2 . If we use the fact 



that F z (zi,z 2 ) 
of the form 



F z (z\ + z 2 ) we can define a polynomial 



F z (z 1 +z 2 ) 



E-9i( z i +z 2 y 

i=0 



(18) 



which contains only P z + 1 unknown coefficients gi. In 
this case we should consider not two separate spectra of 
dynamic variables z\ and z 2 but one united spectrum 
z = z\ + z 2 which is actually measured. As a result, it is 
possible to obtain a tip-surface interaction spectrum F z 



F z = x'^hz 



(19) 



FIG. 4. (Color online) Total transfer function \ (blue solid 
line) for the cantilever with two fiexural modes with trans- 
fer functions xi an d X2 (light and dark green dotted lines 
respectively) . 



are discarded. We also require that the second resonance 
frequency uj 2 is not an integer multiplier of the cj±, so 
the second band does not capture any higher harmon- 
ics and intermodulation products produced by the drive 
in the first band. If ui 2 were a harmonic of lj\ it would 
be of considerable advantage for force measurement!^. 
Considering the monomial basis ( 18 ) we can approxi- 



mately evaluate the Fourier spectrum of the i th power 
J-[z k ] = F[z] * ■ ■ ■ * F[z] via convolution of the transfer- 
function x with itself, which gives an upper bound of the 
response in the frequency domain. Figure [5]ja) demon- 
strates that only components with odd powers i have sig- 
nificant value in the narrow bands near the resonances. 
Consequentially, only parameters gi for odd powers z l (i 
odd) can be found directly using the measured spectrum. 

Having solved the system for the odd parameters, we 
can use them to recover the even parameters by applying 
an additional constraint: that the tip-surface force for tip 
positions above its rest point equals zero F z (z > 0) = 
which leads to the system 



N/ 2 

i=0 i=0 



N/2 



v 2i+l 



(20) 



for all Zfc > measured at discrete time moments t^. 

Reconstruction of the velocity-dependent part of F z is 
achieved by adding a new variable z — z.\ + z 2 to the 
model ( 18 ) which consequently increases the number of 



making use of the total response functional x = Xl + X2 parameters by the degree of the polynomial in i 
depicted in fig. [3] 

Prior to force reconstruction it is necessary to investi- 
gate what kind of information about F z is contained in 
the spectral bands for z\ and z 2 . In the current investi- 
gation we try to reconstruct the force using information 
contained only in the narrow frequency bands near u>i 
and lu 2 , so all weak response peaks outside these bands 



£(*.*) = EE 



gij* 



(21) 



i=0 j=0 



Noticing that T\z l 0\ oc F[z t+ i], only coefficients <?y in 
front of z l z 3 where i+j is odd can be determined from the 
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FIG. 5. (Color online) (a) Response spectra for powers of single dynamic variable with two resonances and (b) two separate 
variables with one resonance. Two frequency windows near resonances wi and 0J2 are highlighted with light and dark green 
colors respectively. Red dashed line means that the corresponding maximum response lies outside the frequency bands. 



measured spectrum. While coefficients before 



zz , etc. can be found using the system (20). 



We simulate the collinear bimodal case using the 
CVODE integrator^ 1 with: wi = 2tt 300 kHz, ki = 40 
N/m, Qi = 400, uj 2 = 6.3wi, k 2 = 50fci^Q 2 = 3Qi (ra- 
tios for the second mode are taken frorrP^l). The driving 
forces f 1 2 are chosen to have the same phase and give 
equal maximum free response (when F z = 0) at each 
mode A Zl — A Z2 = 12.5 nm so the total maximum am- 
plitude of oscillations is A z = 25 nm; all four drive fre- 
quencies W1.2 i Suj/2 are integer multipliers of base fre- 
quency 5u = 27r0.2 kHz. The engaged height h above 
the surface is 17 nm. The model of the tip-surface force 
F z is the vdW-DMT force^ with the nonlinear damping 
term exponentially dependent on the tip position 



z > o,q 



F z {z, z) = F c z on (z) + Ff s {z, z) 

{HR 
-If + tW^o-*), z<a 1221 
Ff s (z,z) = - ( 7l i + 73 i 3 ) e - z / x * 



with the following seven phenomenological parameters: 
intermolecular distance clq = 0.3 nm, Hamaker constant 
H = 7.1 x 10" 20 J, effective modulus E* — 1.0 GPa, tip 
radius R = 10 nm, damping decay length X z — 1.5 nm, 
damping factors 71 = 2.2 x 10~ 7 kg/s and 73 = 10~ 22 
kg-s/m 2 (fig^a)). It is worth noting that the calibration 
of the higher cigenmodes parameters uji , Qi and ki which 
is required for force reconstruction, is itself a challenging 
task in multimodal AF]VpH^l. 

Using B12 — 24 peaks in each band (fig. [3| for the 
reconstruction, we assume the model (21) degree in z 
to be P z = 21. Numerical results show that the force 
reconstruction using higher powers of z l (i > 1) is less 
reliable as it encounters difficulties connected with multi- 
variate interpolation. Although two modes give us twice 
the number of spectral components in comparison with 
singlemode case, this additional information is still in- 
sufficient for reconstruction of nonlinear (non-viscous) 



damping. Therefore, we are restricted the model to be 
linear in i. 

This numerical analysis suggests that there is no dif- 
ference in the quality of the reconstructed force using in- 
termodulation products from frequency bands surround- 
ing the first, second or both resonances. We show re- 
sults for the reconstructed F z (z, z) and its cross-sections 
in fig. j6^b) and fig. [7] using only the resonant detection 
band around the first eigenmode. The linear fit (21) for 



the model with nonlinear damping ( 22 ) nicely captures 



the overall trend of the dissipative part and the recon- 
structed conservative part demonstrates excellent agree- 
ment with the actual force. If we simulate the system 



(16) using the force F z linear in i, for instance assuming 
73 = in (22), the reconstructed force F z shows nearly 



perfect agreement with the actual force F z (figures are 
not included). 

One can compare the information contained in the two 
frequency bands near each eigenmode resonance by esti- 
mating the quality of the reconstruction as a function 
of the number of spectral components B12 used in the 
reconstruction, as the least square error function 



(23) 



{Bi)= / F z (z,z\Bi)-F z {z,z) 



This error function plotted versus the number of spectral 
components is depicted in fig. [8j We see similar behav- 
ior for both bands: significant drop for the number of 
spectral components lager than half of the polynomial 
power P z = 21 in the expansion of F z , and no qualita- 
tive improvements after reaching a number of spectral 
components less than P z . 

Another interesting observation regards the recon- 
struction of a dissipative force which turns-on only for 
the tip velocities i, higher than the maximum velocity 
of the first mode Z\ (for some constant amplitude) and 
is zero otherwise. Information from both bands gives 
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z (nm) z (nm) 

FIG. 6. (Color online) (a) vdW-DMT force with position dependent damping in the region of the free tip motion with bimodal 
drive. White and black dashed lines confine covered domains of phase space in case of single mode drive of the first and second 
modes respectively (difference between maximum velocities is about order of magnitude). In all three cases, total maximum 
amplitude of oscillations is kept constant and equals 25 nm. (b) Reconstructed force in the region of the engaged tip motion. 
Its cross-sections 1-3 are depicted in fig. [TJa)-(c) to highlight agreement between the actual force used in simulation and the 
reconstructed force. 




z (nm) z (nm) z (nm/s) 

FIG. 7. (Color online) Cross-sections of the reconstructed tip-surface force (black) with comparison to the actual force used 
in simulation (red) for the cantilever driven at two flexural modes, (a) F z (z,i — 0); (b) F z (z,z = 0.2i max ) — F z (z,z = 0); 
(c) F z (z = 0.75z m i n ,i). Linear fit for nonlinear damping catches nicely overall trend and conservative part demonstrates an 
excellent agreement with the actual force. 




number of spectral components number of spectral components 

FIG. 8. (Color online) Absolute error of reconstruction of (a) conservative part of the tip-surface force (fig. |7[a)) and (b) 
dissipative part (fig. [Fpa)) versus number of spectral components taken into account in the first frequency band (blue circles) 
and in the second (green squares). 



approximately the same reconstructed curves with show- 
ing overall trend of the dissipative part and in excellent 
agreement for the conservative part (plots not shown). 
However, reconstructing with spectral components from 
only the first frequency band yields more accurate dissi- 
pative force approximation as Z\ has smaller stored oscil- 
lation energy and, consequently, is more vulnerable to the 
dissipative force than Z2- Nonetheless, if we excite only 
the first mode, we would not be able to reconstruct this 
"threshold" dissipative force at all, as the magnitude of 
the tip velocity would not be enough to turn on the dissi- 
pation. Thus, simultaneous excitation of two eigenmodes 
allows to explore a wider region of the phase space of the 
tip motion (fig. [6^ a)) while keeping the total maximum 
amplitude constant. 



B. Case 2: Two orthogonal modes 

This case corresponds to the dynamics of one flexural 
z and one torsional y mode 



k * + Q7^ i + Z ) =F z (z,y)+f z (t) 

k v {^y + o^y + y)= F v& v) + 



(24) 



with two different projections of the conservative force, 
F z and F y on the r.h.s. 

We start from the polynomial models for reconstruc- 
tion 



Pz Py 



F z (z,y) = E Eff£ } *V 

i=0 j=0 
i=0j=0 



(25) 



As for case 1 it is not possible to find all parameters 
of these models. Firstly, as we are limited in number 
of measurable intermodulation products and therefore in 
maximum degree of the polynomial. We choose z direc- 
tion as the most interesting degree of freedom, by which 
we mean that the maximum degree of the polynomial in 
this variable will be much higher than for y. In accor- 
dance to the fig. [5jb), the captured information about 
forces F z will be odd in z and even in y and vice versa 



for F, 



u 



F z (-z,±y) = -F z (z,y) 
F y (±z,-y) = 



F y {z,y) 



(26) 



as the first flexural resonance uj z is typically far lower 
in frequency than the first torsional resonance tuj^. It is 

(z) 

possible to recover the coefficients of even powers of y 

and gf^ of odd powers of y (when i+j is even) by using 
the additional constraint for z dependence of the force 
components F z , y (z > 0,y) — and eq. (20). While the 



of y and F y with even powers of y is lost because we have 
no such constraint on the y dependence. 

Simulation parameters for the system (24 1 are: uj z = 
2tt300 kHz, k z = 40 N/m, Q z = 400, uj y = 6.3uj y , k y = 
50k y , Qi = 3Q y . The driving forces f Zj!/ are chosen to 
have the same phase and give maximum free response 
(when F z y = 0) A z — 25 nm and A y = 12.5 nm; all four 
drive frequencies u> z ^ y ± Suj/2 are integer multipliers of 
base frequency Slj = 27r0.2 kHz. The engaged height h 
above the surface is 17 nm. The model for the component 
of the tip-surface force perpendicular to the surface is the 
same vdW-DMT force (22) used in the previous case, 
without the dissipation term F* s . The model of the 
force component parallel to the surface is a nonlinear 
conservative restoring force 



F y (z,y) = -(c iy + c 3 y 3 ) e'*'^ 



(27) 



where X z = 1.5 nm, c\ = 0.22 N/m and C3 = 0.1 N/m 3 
are constants. These two components of F ts (r) are illus- 
trated in fig. [9]ja)-(b). 

Using only 24 intermodulation peaks in each band for 
z and y, the spectral fitting method reconstructs the two- 
dimensional vector force field F ts ( 22 ) and ([27]) up to the 
21 st power in z and third power in y (fig. I9fc)-(d) and 
10). The reconstructed force is in good agreement with 



information about all coefficients of F z with odd powers 



the actual model and perfect agreement is reachable if we 
assume the model for F z as eq. (25) only, independent 
of y. 



V. CONCLUDING REMARKS 

In this paper we have discussed the basic problem 
of multimodal ImAFM and we proposed a theoretical 
framework for reconstructing multidimensional forces us- 
ing this technique. We demonstrated the possibility of 
reconstructing tip-surface interactions for two character- 
istic bimodal cases. We have studied tip-surface interac- 
tions that show the most nonlinear behavior in the first 
degree of freedom, and are linear or cubic in the sec- 
ond degree of freedom. As the numerical results have 
shown, it is possible to accurately reconstruct dependen- 
cies up to a cubic order in the second degree of freedom 
for the two-dimensional force model using information 
from only one frequency band. We found that excitation 
of two flexural modes with two well separated resonances 
does not allow for a precise reconstruction of a nonlinear 
damping force using only the information contained in 
corresponding narrow bands near the resonances. How- 
ever, the reconstruction nicely captures the overall linear 
trend, or that of a viscous damping. The use of a second, 
higher frequency eigenmode allows for reconstruction on 
a wider region of the phase space of a tip motion, en- 
abling exploration of dissipative interactions inaccessible 
to the first mode alone for a given maximum amplitude 
of motion. Additionally, the first eigenmode is found to 
be more sensitive to dissipation forces acting on the tip. 
Finally, using simultaneous excitation of two orthogonal 



9 




-25 -20 -15 -10 



F(nN) 



(b) 



§ 



-25 -20 -15 -10 



z (nm) 



z (nm) 




50.00 
25.00 
0.000 



™ -50.00 

F(nN) 




H 50.00 
25.00 
0.000 
_ -25.00 
^ -50.00 

F(nN) 



FIG. 9. (Color online) Components F z and .F y of two-dimensional conservative force in the region of the free tip motion (a,b) 
and reconstructed components in the region of the engaged tip motion (c,d). Cross-sections 1-3 are illustrated in fig. 10 a)-(c) 
to highlight agreement between the actual force used in simulation and the reconstructed force. 
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FIG. 10. (Color online) Cross-sections of reconstructed tip-surface force (black) with comparison to the actual force used in 
simulation (red) for the cantilever driven at two orthogonal modes, (a) F z (z); (b) F y (z,y = 0.2i max ); (c) F y (z = 0.75z m i n , y). 
Reconstructed force is in a good agreement with underlying model. Perfect agreement is reachable if we assume model for F z 
independent of y. 



modes we can reconstruct nonlinear position-dependent 
lateral forces simultaneously with vertical forces. This 
approach represents a path toward the determination of 
a vectorial force field by frequency domain multiplexing 
of the multimodal response of an AFM cantilever. 



Appendix: Generalized eigencoordinates 



We start from the governing equation for a two- 
dimensional cantilever 



(Gxy + St) [w (x, y, t)] = F(x, y, t) 



(A.l) 
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and try to find solution w (x, y, t) separated in time and 
space, and expanded into the set of normal modes 



are effective stiffness and mass of corresponding degree 
of freedom, and considering damping and inertia 



(A.2) 



i=0 



Functions <pi form orthonormal set on the geometrical 
shape of cantilever J7 C 



(A.3) 



where 8j is Kronecker delta equals 1 for i = j and 
otherwise. 



Inserting solution (A.2) into eq. (A.l) with following 



multiplication by <j>i and integrating over the plane f2 c 
yields a system of differential equations for the general- 
ized coordinates qt(t) 



Qi J <t>iG X y(t>idfL X y+gt[qi] J 4>i dCl xy — J F{Q,,t)(j>ij<iQ. xy 

(A.4) 



here the orthonormal condition (A.3) is used 



Denoting Qi as a differential operator governing motion 
of i th generalized coordinate 



where 



Qi = ki + rriiQ t 



h = J 4>iQ X y4>i d^a 

nii = J $ dQ xy 



(A.5) 



(A.6) 



d 2 



d 



dt 2 + 1 ~dt 



(A.7) 



with some constant 7 (homogeneous viscous medium 
damping), we arrive at the final system 



1 .. 1 



M , ,2 



2 ^* 1 /"i 



qi + ~ — qi +qi) =Fi (t) + U (t) (A. 



where resonant frequencies wj — yjki/rrii and quality fac- 
tors Qi — uji /j are introduced. Here the time-dependent 
forces 



F, 



(t) + U (t) := J faix, y)F(x, y, t) dfl xy (A.9) 



represent anharmonic contribution of the tip-surface in- 
teraction and harmonic contribution of the drive, respec- 
tively. 
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